rm(list = ls(all = TRUE))

gc()
gc()

set.seed(123456)

1 libraries

library(scales)
library(ggplot2)
# packageVersion("ggplot2")
library(magrittr)
# library(crayon)
library(rstatix)
# library(ggpubr)
# library(Compositional)
# library(equalCovs)
# library(psych)
library(corrplot) 
# library(pdist)
library(RColorBrewer)
library(plotly)
library(heatmaply)
# library(Hmisc)
library(htmlwidgets)
library(viridis)
library(ggh4x)


`%nin%` = Negate(`%in%`)

2 f-ctions

3 predefined contrasts

fpi = file.path('..', 'input')
fn = 'contrasts.txt'

comp = data.table::fread(file.path(fpi, fn), header = FALSE)

4 metabolomics

dir.create(file.path('..', 'reports', 'metabolomics'))

cnt = 1
fpi = file.path('..', 'input')


fn = "data_targeted.xlsx"
omics = openxlsx::read.xlsx(file.path(fpi, fn), sheet = 'metabolomics-leaves')
data.table::setDF(omics)

Genotype = 'Desiree'
omics = cbind(Genotype, omics)


PostInductionDay = as.numeric(stringr::str_split_i(gsub('S', '', omics$SampleID), '_', i = 2))
omics = cbind(PostInductionDay, omics)

omics = omics[omics$PostInductionDay %in% 1:14, ]

Condition = gsub('_.*', '', omics$SampleID)
omics = cbind(Condition, omics)
omics = omics[-which(omics$Condition == 'W' & omics$PostInductionDay %in% c(8,14)), ]

omics$Condition = factor(omics$Condition, levels = c('C', 'H', 'D', 'HD', 'W'))

table(omics$Condition, omics$PostInductionDay)
##     
##      1 7 8 14
##   C  4 4 4  4
##   H  4 4 4  4
##   D  0 0 4  4
##   HD 0 0 4  4
##   W  4 4 0  0
mygroup = paste(omics$Condition, stringr::str_pad(omics$PostInductionDay, 1, pad = "0"), sep = '_')
# table(mygroup %in% comp$V1)
# table(mygroup %in% comp$V3)

omics = cbind(mygroup, omics)

4.1 palette

n = 6


# my.ggplot.palette = brewer.pal(8, "Set2")[c(8, 2, 6, 4, 5, 3)]#[c(8, 2, 6, 7, 5, 3)]
my.ggplot.palette = c('#7F7F7F', '#C00000', '#FFC000', '#ED7D31',  
                      # '#70AD47', 
                      '#4472C4')
pie(rep(1, length(my.ggplot.palette)), 
    col = my.ggplot.palette, 
    main = 'ggplot palette',
    labels = c('C', 'H', 'D', 'HD', 
               # 'HDW', 
               'W'))

extend.palette = cbind(my.ggplot.palette, c('C', 'H', 'D', 'HD', 
                                            # 'HDW', 
                                            'W'))
colnames(extend.palette) = c('colour', 'group')
extend.palette = merge(extend.palette, omics[, 1:4], by.x = 'group', by.y = 'Condition', all = TRUE)


par(mfrow = c(2,2))

nn = 9
# (intervals = cut(c(-1,1), n*2))
intervals = cut(c(-1,1), 
                breaks = setdiff(seq(-1, 1, 0.1), 0), 
                include.lowest = TRUE)

my.heatmap.col1 = c(rev(brewer.pal(nn, 'Reds')), 'white', brewer.pal(nn, 'Blues'))
# my.heatmap.col1 = cividis(19, direction = -1)
pie(rep(1, length(my.heatmap.col1)), 
    col = rev(my.heatmap.col1), 
    labels = levels(intervals), 
    main = 'heatmap palette 1')


my.heatmap.col2 = c(my.heatmap.col1[1:3], 
                    rep('grey90', length(4:(nn*2-2))),
                    my.heatmap.col1[(nn*2-1): (nn*2 + 1)])
pie(rep(1, length(my.heatmap.col1)), 
    col = rev(my.heatmap.col2), 
    labels = levels(intervals), 
    main = 'heatmap palette 2')


par(mfrow = c(1,1))

4.2 density

# column in which measurements start
nx = 6

data = omics
# colnames(data)[nx:ncol(data)]

data[, nx:ncol(data)] = apply(data[, nx:ncol(data)] , 2, as.numeric)

# https://www.cookbook-r.com/Manipulating_data/Converting_data_between_wide_and_long_format/
data.long = tidyr::gather(data, molecule, measurement, colnames(data)[nx]:colnames(data)[ncol(data)], factor_key=FALSE)


data.long$measurement = as.numeric(data.long$measurement)

data.long$molecule = factor(data.long$molecule, levels = colnames(data)[nx:ncol(data)])

data.long$PostInductionDay = ordered(data.long$PostInductionDay, levels = sort(unique(data.long$PostInductionDay )))



data.long = data.long[!is.na(data.long$measurement), ]


# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = Condition, colour = Condition)) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = my.ggplot.palette) +
  scale_fill_manual(name = 'Legend',
                    values = my.ggplot.palette) 

# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = mygroup , colour = mygroup )) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = unique(extend.palette[, 2:3])[, 1]) +
  scale_fill_manual(name = 'Legend',
                    values = unique(extend.palette[, 2:3])[, 1]
                    ) 

data.long.SE = summarySE(data = data.long, 
                         measurevar="measurement", 
                         groupvars=c("molecule", "Genotype", "Condition", "PostInductionDay", "mygroup"), 
                         na.rm = TRUE)
ncol = 6

jitter = position_jitter(width = 0.0, height = 0.0)
dodge = position_dodge(width=0.0)

4.3 plot measurements

cd = sort(unique(data.long.SE$PostInductionDay))

mypallette = my.ggplot.palette

# useful in case of more genotypes
lims = data.long.SE %>% 
  dplyr::group_by(molecule) %>% 
  dplyr::summarise(Panel = molecule, 
                   ymin=min(mean-se),ymax=max(mean+se),
                   n = 5
                   )
lims = unique(lims)
lims = split(lims, lims$Panel) # make a list
lims = lapply(lims, function(l) {
  scale_y_continuous(limits = c(l$ymin, l$ymax),
                     n.breaks = l$n)
  }
)


title = 'Desiree'
data.SE = data.long.SE[data.long.SE$Genotype == title, ]
myggplot(data.SE, title, cnt, lncol = 2, lims, subdir = 'metabolomics')

4.4 correlation

mytitle = 'Desiree'
df = omics[omics$Genotype == mytitle, ]
df = as.matrix(df[, nx:ncol(df)])
df = apply(df, 2, as.numeric)
cs = colSums(is.na(df))/nrow(df)
if (any(cs == 1)) df = df[, -which(cs == 1)]

rownames(df) = omics$SampleID[omics$Genotype == mytitle]

mycond = levels(omics$Condition)

# X11();
# plot.new();
my.pairs.panels(df, mytitle, mycond, cccex = 1)

dev.copy(pdf,
         file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   3
# # try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()


my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')





mycond = 'C'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# # try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')


mycond = 'H'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')


mycond = 'D'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')


mycond = 'W'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')

    
mycond = 'HD'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'metabolomics')
genotype = 'Desiree'
De.log2FC = NULL
De.p = NULL

for (i in 1:nrow(comp)) {
  
  # print(i)
  # print( as.character(comp[i, ]))
  
  dtX = data.long[data.long$mygroup %in% comp[i, c(1,3)] & data.long$Genotype %in% genotype, ]

  stress = dtX[dtX$mygroup == as.character(comp[i, 1]), ]
  Ctrl = dtX[dtX$mygroup == as.character(comp[i, 3]), ]
  
  StressC = gsub('_.*', '', comp[i, 1])
  CtrlC = gsub('_.*', '', comp[i, 3])
  
  log2FC = my.logFC(Ctrl = Ctrl, stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  p = my.customised.t.test(Ctrl = Ctrl, Stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  
  De.log2FC = rbind(De.log2FC, log2FC)
  De.p = rbind(De.p, p)
  
  
}


Desiree.stat = merge(De.log2FC, De.p,
                     by.y = c("Genotype", "molecule", "group2", "group1"),
                     by.x = c("Genotype", "molecule", "stressGroup", "controlGroup"),
                     all = TRUE)


openxlsx::write.xlsx(Desiree.stat, '../reports/Desiree.logFC-metabolomics_leaves.xlsx', asTable = TRUE)

4.5 heatmap

comp = as.data.frame(comp)
forHeatmap = comp[intersect(grep('C', comp[, 1], invert = TRUE), grep('C', comp[, 3])), ]
forHeatmap$comparison = paste(forHeatmap[, 1], forHeatmap[, 3], sep = '-')

df = Desiree.stat[Desiree.stat$stressGroup %in% forHeatmap[, 1] & Desiree.stat$controlGroup %in% forHeatmap[, 3], ]

my.pheatmap.logFC(df, forHeatmap,
                  main = 'Desiree-metabolomics',
                  cluster_rows = FALSE, cluster_cols = FALSE,
                  subdir = 'metabolomics',
                  width = 6, height = 6, cnt = cnt)
# rm()
gc()
gc()

5 hormonomics

dir.create(file.path('..', 'reports', 'hormonomics'))

cnt = 2
fpi = file.path('..', 'input')


fn = "data_targeted.xlsx"
omics = openxlsx::read.xlsx(file.path(fpi, fn), sheet = 'hormonomics-leaves')
data.table::setDF(omics)

Genotype = 'Desiree'
omics = cbind(Genotype, omics)


PostInductionDay = as.numeric(stringr::str_split_i(gsub('S|P', '', omics$SampleID), '_', i = 2))
omics = cbind(PostInductionDay, omics)

omics = omics[omics$PostInductionDay %in% 1:14, ]

Condition = gsub('_.*', '', omics$SampleID)
omics = cbind(Condition, omics)
omics = omics[-which(omics$Condition == 'W' & omics$PostInductionDay %in% c(8,14)), ]

omics$Condition = factor(omics$Condition, levels = c('C', 'H', 'D', 'HD', 'W'))

table(omics$Condition, omics$PostInductionDay)
##     
##      1 7 8 14
##   C  4 4 4  4
##   H  4 4 4  4
##   D  0 0 4  4
##   HD 0 0 4  4
##   W  4 4 0  0
mygroup = paste(omics$Condition, stringr::str_pad(omics$PostInductionDay, 1, pad = "0"), sep = '_')
# table(mygroup %in% comp$V1)
# table(mygroup %in% comp$V3)

omics = cbind(mygroup, omics)

5.1 palette

my.ggplot.palette = c('#7F7F7F', '#C00000', '#FFC000', '#ED7D31',  
                      # '#70AD47', 
                      '#4472C4')

extend.palette = cbind(my.ggplot.palette, c('C', 'H', 'D', 'HD',
                                            #'HDW', 
                                            'W'))
colnames(extend.palette) = c('colour', 'group')
extend.palette = merge(extend.palette, omics[, 1:4], by.x = 'group', by.y = 'Condition', all = TRUE)

5.2 density

# column in which measurements start
nx = 6

data = omics
# colnames(data)[nx:ncol(data)]

data[, nx:ncol(data)] = apply(data[, nx:ncol(data)] , 2, as.numeric)

data.long = tidyr::gather(data, molecule, measurement, colnames(data)[nx]:colnames(data)[ncol(data)], factor_key=FALSE)


data.long$measurement = as.numeric(data.long$measurement)

data.long$molecule = factor(data.long$molecule, levels = colnames(data)[nx:ncol(data)])

data.long$PostInductionDay = ordered(data.long$PostInductionDay, levels = sort(unique(data.long$PostInductionDay )))


data.long = data.long[!is.na(data.long$measurement), ]


# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = Condition, colour = Condition)) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = my.ggplot.palette) +
  scale_fill_manual(name = 'Legend',
                    values = my.ggplot.palette) 

# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = mygroup , colour = mygroup )) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = unique(extend.palette[, 2:3])[, 1]) +
  scale_fill_manual(name = 'Legend',
                    values = unique(extend.palette[, 2:3])[, 1]
                    ) 

data.long.SE = summarySE(data = data.long, 
                         measurevar="measurement", 
                         groupvars=c("molecule", "Genotype", "Condition", "PostInductionDay", "mygroup"), 
                         na.rm = TRUE)

5.3 plot measurements

cd = sort(unique(data.long.SE$PostInductionDay))

mypallette = my.ggplot.palette


lims = data.long.SE %>% 
  dplyr::group_by(molecule) %>% 
  dplyr::summarise(Panel = molecule, 
                   ymin=min(mean-se),ymax=max(mean+se),
                   n = 5
                   )
lims = unique(lims)
lims = split(lims, lims$Panel) # make a list
lims = lapply(lims, function(l) {
  scale_y_continuous(limits = c(l$ymin, l$ymax),
                     n.breaks = l$n)
  }
)


title = 'Desiree'
data.SE = data.long.SE[data.long.SE$Genotype == title, ]
myggplot(data.SE, title, cnt, lncol = 2, lims, subdir = 'hormonomics')

omics  = omics[, -grep('neoPA', colnames(omics))]
data.long.SE  = data.long.SE[data.long.SE$molecule != 'neoPA', ]
data.long  = data.long[data.long$molecule != 'neoPA', ]

5.4 correlation

mytitle = 'Desiree'
df = omics[omics$Genotype == mytitle, ]
df = as.matrix(df[, nx:ncol(df)])
df = apply(df, 2, as.numeric)
cs = colSums(is.na(df))/nrow(df)
if (any(cs == 1)) df = df[, -which(cs == 1)]

rownames(df) = omics$SampleID[omics$Genotype == mytitle]

mycond = levels(omics$Condition)

# X11();
# plot.new();
my.pairs.panels(df, mytitle, mycond, cccex = 1)

dev.copy(pdf,
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##  11
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()


my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')





mycond = 'C'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')


mycond = 'H'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')


mycond = 'D'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')


mycond = 'W'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')

    
mycond = 'HD'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'hormonomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'hormonomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'hormonomics')
genotype = 'Desiree'
De.log2FC = NULL
De.p = NULL

for (i in 1:nrow(comp)) {
  
  # print(i)
  # print( as.character(comp[i, ]))
  
  dtX = data.long[data.long$mygroup %in% comp[i, c(1,3)] & data.long$Genotype %in% genotype, ]

  stress = dtX[dtX$mygroup == as.character(comp[i, 1]), ]
  Ctrl = dtX[dtX$mygroup == as.character(comp[i, 3]), ]
  
  StressC = gsub('_.*', '', comp[i, 1])
  CtrlC = gsub('_.*', '', comp[i, 3])
  
  log2FC = my.logFC(Ctrl = Ctrl, stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  p = my.customised.t.test(Ctrl = Ctrl, Stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  
  De.log2FC = rbind(De.log2FC, log2FC)
  De.p = rbind(De.p, p)
  
  
}


Desiree.stat = merge(De.log2FC, De.p,
                     by.y = c("Genotype", "molecule", "group2", "group1"),
                     by.x = c("Genotype", "molecule", "stressGroup", "controlGroup"),
                     all = TRUE)


openxlsx::write.xlsx(Desiree.stat, '../reports/Desiree.logFC-hormonomics_leaves.xlsx', asTable = TRUE)

5.5 heatmap

comp = as.data.frame(comp)
forHeatmap = comp[intersect(grep('C', comp[, 1], invert = TRUE), grep('C', comp[, 3])), ]
forHeatmap$comparison = paste(forHeatmap[, 1], forHeatmap[, 3], sep = '-')

df = Desiree.stat[Desiree.stat$stressGroup %in% forHeatmap[, 1] & Desiree.stat$controlGroup %in% forHeatmap[, 3], ]

my.pheatmap.logFC(df, forHeatmap,
                  main = 'Desiree-hormonomics',
                  cluster_rows = FALSE, cluster_cols = FALSE,
                  subdir = 'hormonomics',
                  width = 6, height = 6, cnt = cnt)
# rm()
gc()
gc()

6 transcriptomics

dir.create(file.path('..', 'reports', 'transcriptomics'))

cnt = 3
fpi = file.path('..', 'input')


fn = "data_targeted.xlsx"
omics = openxlsx::read.xlsx(file.path(fpi, fn), sheet = 'transcriptomics-leaves')
data.table::setDF(omics)

Genotype = 'Desiree'
omics = cbind(Genotype, omics)


PostInductionDay = as.numeric(stringr::str_split_i(gsub('S|P', '', omics$SampleID), '_', i = 2))
omics = cbind(PostInductionDay, omics)

omics = omics[omics$PostInductionDay %in% 1:14, ]

Condition = gsub('_.*', '', omics$SampleID)
omics = cbind(Condition, omics)
omics = omics[-which(omics$Condition == 'W' & omics$PostInductionDay %in% c(8,14)), ]

omics$Condition = factor(omics$Condition, levels = c('C', 'H', 'D', 'HD', 'W'))

table(omics$Condition, omics$PostInductionDay)
##     
##      1 7 8 14
##   C  4 4 4  4
##   H  4 4 4  4
##   D  0 0 4  4
##   HD 0 0 4  4
##   W  4 4 0  0
mygroup = paste(omics$Condition, stringr::str_pad(omics$PostInductionDay, 1, pad = "0"), sep = '_')
# table(mygroup %in% comp$V1)
# table(mygroup %in% comp$V3)

omics = cbind(mygroup, omics)

6.1 palette

my.ggplot.palette = c('#7F7F7F', '#C00000', '#FFC000', '#ED7D31',  
                      # '#70AD47', 
                      '#4472C4')

extend.palette = cbind(my.ggplot.palette, c('C', 'H', 'D', 'HD',
                                            #'HDW', 
                                            'W'))
colnames(extend.palette) = c('colour', 'group')
extend.palette = merge(extend.palette, omics[, 1:4], by.x = 'group', by.y = 'Condition', all = TRUE)

6.2 density

# column in which measurements start
nx = 6

data = omics
# colnames(data)[nx:ncol(data)]

data[, nx:ncol(data)] = apply(data[, nx:ncol(data)] , 2, as.numeric)

data.long = tidyr::gather(data, molecule, measurement, colnames(data)[nx]:colnames(data)[ncol(data)], factor_key=FALSE)


data.long$measurement = as.numeric(data.long$measurement)

data.long$molecule = factor(data.long$molecule, levels = colnames(data)[nx:ncol(data)])

data.long$PostInductionDay = ordered(data.long$PostInductionDay, levels = sort(unique(data.long$PostInductionDay )))


data.long = data.long[!is.na(data.long$measurement), ]


# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = Condition, colour = Condition)) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = my.ggplot.palette) +
  scale_fill_manual(name = 'Legend',
                    values = my.ggplot.palette) 

# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = mygroup , colour = mygroup )) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = unique(extend.palette[, 2:3])[, 1]) +
  scale_fill_manual(name = 'Legend',
                    values = unique(extend.palette[, 2:3])[, 1]
                    ) 

data.long.SE = summarySE(data = data.long, 
                         measurevar="measurement", 
                         groupvars=c("molecule", "Genotype", "Condition", "PostInductionDay", "mygroup"), 
                         na.rm = TRUE)

6.3 plot measurements

cd = sort(unique(data.long.SE$PostInductionDay))

mypallette = my.ggplot.palette


lims = data.long.SE %>% 
  dplyr::group_by(molecule) %>% 
  dplyr::summarise(Panel = molecule, 
                   ymin=min(mean-se),ymax=max(mean+se),
                   n = 5
                   )
lims = unique(lims)
lims = split(lims, lims$Panel) # make a list
lims = lapply(lims, function(l) {
  scale_y_continuous(limits = c(l$ymin, l$ymax),
                     n.breaks = l$n)
  }
)


title = 'Desiree'
data.SE = data.long.SE[data.long.SE$Genotype == title, ]
myggplot(data.SE, title, cnt, lncol = 2, lims, subdir = 'transcriptomics')

6.4 correlation

mytitle = 'Desiree'
df = omics[omics$Genotype == mytitle, ]
df = as.matrix(df[, nx:ncol(df)])
df = apply(df, 2, as.numeric)
cs = colSums(is.na(df))/nrow(df)
if (any(cs == 1)) df = df[, -which(cs == 1)]

rownames(df) = omics$SampleID[omics$Genotype == mytitle]

mycond = levels(omics$Condition)

# X11();
# plot.new();
my.pairs.panels(df, mytitle, mycond, cccex = 1)

dev.copy(pdf,
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##  19
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()


my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')





mycond = 'C'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')


mycond = 'H'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')


mycond = 'D'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new();
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')


mycond = 'W'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')

    
mycond = 'HD'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'transcriptomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'transcriptomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'transcriptomics')
genotype = 'Desiree'
De.log2FC = NULL
De.p = NULL

for (i in 1:nrow(comp)) {
  
  # print(i)
  # print( as.character(comp[i, ]))
  
  dtX = data.long[data.long$mygroup %in% comp[i, c(1,3)] & data.long$Genotype %in% genotype, ]

  stress = dtX[dtX$mygroup == as.character(comp[i, 1]), ]
  Ctrl = dtX[dtX$mygroup == as.character(comp[i, 3]), ]
  
  StressC = gsub('_.*', '', comp[i, 1])
  CtrlC = gsub('_.*', '', comp[i, 3])
  
  log2FC = my.logFC(Ctrl = Ctrl, stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  p = my.customised.t.test(Ctrl = Ctrl, Stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  
  De.log2FC = rbind(De.log2FC, log2FC)
  De.p = rbind(De.p, p)
  
  
}


Desiree.stat = merge(De.log2FC, De.p,
                     by.y = c("Genotype", "molecule", "group2", "group1"),
                     by.x = c("Genotype", "molecule", "stressGroup", "controlGroup"),
                     all = TRUE)


openxlsx::write.xlsx(Desiree.stat, '../reports/Desiree.logFC-transcriptomics_leaves.xlsx', asTable = TRUE)

6.5 heatmap

comp = as.data.frame(comp)
forHeatmap = comp[intersect(grep('C', comp[, 1], invert = TRUE), grep('C', comp[, 3])), ]
forHeatmap$comparison = paste(forHeatmap[, 1], forHeatmap[, 3], sep = '-')

df = Desiree.stat[Desiree.stat$stressGroup %in% forHeatmap[, 1] & Desiree.stat$controlGroup %in% forHeatmap[, 3], ]

my.pheatmap.logFC(df, forHeatmap,
                  main = 'Desiree-transcriptomics',
                  cluster_rows = FALSE, cluster_cols = FALSE,
                  subdir = 'transcriptomics',
                  width = 6, height = 6, cnt = cnt)
# rm()
gc()
gc()

7 tubers metabolomics

dir.create(file.path('..', 'reports', 'tubers-metabolomics'))

cnt = 4
fpi = file.path('..', 'input')


fn = "data_targeted.xlsx"
omics = openxlsx::read.xlsx(file.path(fpi, fn), sheet = 'metabolomics-tubers')
data.table::setDF(omics)

Genotype = 'Desiree'
omics = cbind(Genotype, omics)


PostInductionDay = as.numeric(stringr::str_split_i(gsub('S|P', '', omics$SampleID), '_', i = 2))
omics = cbind(PostInductionDay, omics)


Condition = gsub('_.*', '', omics$SampleID)
omics = cbind(Condition, omics)

omics$Condition = factor(omics$Condition, levels = c('C', 'H', 'D', 'HD', 'W'))

table(omics$Condition, omics$PostInductionDay)
##     
##      28
##   C  12
##   H  12
##   D  12
##   HD 12
##   W   8
mygroup = paste(omics$Condition, stringr::str_pad(omics$PostInductionDay, 1, pad = "0"), sep = '_')
# table(mygroup %in% comp$V1)
# table(mygroup %in% comp$V3)

omics = cbind(mygroup, omics)

7.1 palette

my.ggplot.palette = c('#7F7F7F', '#C00000', '#FFC000', '#ED7D31',  
                      # '#70AD47', 
                      '#4472C4')

extend.palette = cbind(my.ggplot.palette, c('C', 'H', 'D', 'HD',
                                            #'HDW', 
                                            'W'))
colnames(extend.palette) = c('colour', 'group')
extend.palette = merge(extend.palette, omics[, 1:4], by.x = 'group', by.y = 'Condition', all = TRUE)

7.2 density

# column in which measurements start
nx = 9

data = omics
# colnames(data)[nx:ncol(data)]

data[, nx:ncol(data)] = apply(data[, nx:ncol(data)] , 2, as.numeric)

data.long = tidyr::gather(data, molecule, measurement, colnames(data)[nx]:colnames(data)[ncol(data)], factor_key=FALSE)


data.long$measurement = as.numeric(data.long$measurement)

data.long$molecule = factor(data.long$molecule, levels = colnames(data)[nx:ncol(data)])

data.long$PostInductionDay = ordered(data.long$PostInductionDay, levels = sort(unique(data.long$PostInductionDay )))


data.long = data.long[!is.na(data.long$measurement), ]


# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = Condition, colour = Condition)) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = my.ggplot.palette) +
  scale_fill_manual(name = 'Legend',
                    values = my.ggplot.palette) 

# X11()
ggplot(data.long[data.long$Genotype == 'Desiree', ], aes(log2(measurement), fill = mygroup , colour = mygroup )) +
  geom_density(alpha = 0.1) + 
  ggtitle('Desiree')  +
  theme_classic() +
  scale_colour_manual(name = 'Legend',
                      values = unique(extend.palette[, 2:3])[, 1]) +
  scale_fill_manual(name = 'Legend',
                    values = unique(extend.palette[, 2:3])[, 1]
                    ) 

data.long.SE = summarySE(data = data.long, 
                         measurevar="measurement", 
                         groupvars=c("molecule", "Genotype", "Condition", "PostInductionDay", "mygroup"), 
                         na.rm = TRUE)

7.3 plot measurements

cd = sort(unique(data.long.SE$PostInductionDay))

mypallette = my.ggplot.palette




title = 'Tubers'
lncol = 2
subdir = 'tubers-metabolomics'

# X11();
myplotY = ggplot(data.long, aes(x = Condition , 
                                y = measurement, 
                                group = Condition,
                                shape = Condition,
                                fill = Condition,
                                colour = Condition)) + 
  geom_boxplot() + 
      facet_wrap(~ molecule, ncol = ncol,
             strip.position = "top",
             scales = 'free_y',
             drop = FALSE
             ) +
      ggtitle(title)  +
    theme_classic() +
      scale_colour_manual(name = 'Legend',
                          values = mypallette) +
      scale_fill_manual(name = 'Legend',
                        values = mypallette
                        ) + 
        labs(x = '', 
             y = 'Relative molecule abundance [+/- SE]') + 
        scale_y_continuous(# trans = log2_trans(), 
                           labels = scales::comma
                           ) +

      theme(axis.text.x = element_text (angle = 0, # 45, 
                                        vjust = 1, 
                                        hjust = 1, 
                                        size = 12),
            axis.text.y = element_text (angle = 0, 
                                        vjust = 1, 
                                        hjust = 1, 
                                        size = 12),
            axis.title.x = element_text (angle = 0, # 45, 
                                        #vjust = 1, 
                                        #hjust = 1, 
                                        size = 12,
                                        face = "bold"),
            axis.title.y = element_text (angle = 90, 
                                        # vjust = 1, 
                                        # hjust = 1, 
                                        size = 12,
                                        face = "bold"),
            axis.text = element_text( angle = 90 ),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            panel.border = element_blank(),
            panel.background = element_blank(),
            legend.position = "bottom",
            legend.background = element_rect(fill=NULL, 
                                             size=0.1, linetype="solid"),
            legend.title = element_text(color = NULL,
                                        size = 12,
                                        face = "bold"),
            legend.text = element_text( size = 12),
            strip.text = element_text(size = 12)) + 
  theme(
        legend.box = "horizontal",
        legend.key.width = unit(2, 'cm'),
        legend.position = c(1, 0),
        legend.justification = c(1, 0)) +
  guides(colour = guide_legend(ncol = lncol)) +
  guides(linetype = guide_legend(ncol = lncol)) +
  guides(shape = "none") +
  theme(panel.spacing.y = unit(1.5, "lines")) +
  labs(color  = "Legend", linetype = "Legend", shape = "Legend")

plot(myplotY)  

# dev.off()
  
  
  fng = paste0(stringr::str_pad(cnt, 3, pad = "0"), '_', title, "_L.pdf")
  ggsave(file.path('..', 'reports', subdir, fng), 
         width = 16, height = 12, units = "in")
  fng = paste0(stringr::str_pad(cnt, 3, pad = "0"), '_', title, "_L.svg")
  ggsave(file.path('..', 'reports', subdir, fng),
         width = 16, height = 12, units = "in")

7.4 correlation

mytitle = 'Desiree'
df = omics[omics$Genotype == mytitle, ]
df = as.matrix(df[, nx:ncol(df)])
df = apply(df, 2, as.numeric)
cs = colSums(is.na(df))/nrow(df)
if (any(cs == 1)) df = df[, -which(cs == 1)]

rownames(df) = omics$SampleID[omics$Genotype == mytitle]

mycond = levels(omics$Condition)

# X11();
# plot.new();
my.pairs.panels(df, mytitle, mycond, cccex = 1)

dev.copy(pdf,
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##  27
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()


my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')





mycond = 'C'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')


mycond = 'H'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')


mycond = 'D'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')


mycond = 'W'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')

    
mycond = 'HD'    
df = omics[omics$Genotype == mytitle & omics$Condition == mycond, ]
df = as.matrix(df[, nx:ncol(df)])

# X11(); plot.new(); 
my.pairs.panels(df, mytitle, mycond, cccex = 1)
dev.copy(pdf,
         file.path('..', 'reports', 'tubers-metabolomics', paste0('SPLOM_', 'PCC_', mytitle, '_', paste(mycond, collapse = '-'), ".pdf")),
         width = 13, height = 13)
## pdf 
##   2
# try(dev.off(dev.list()["RStudioGD"]),silent=TRUE); try(dev.off(),silent=TRUE);graphics.off()

my.heatmaply(df, mytitle, mycond, subdir = 'tubers-metabolomics') 
my.pheatmap(df, mytitle, mycond, subdir = 'tubers-metabolomics')
genotype = 'Desiree'
De.log2FC = NULL
De.p = NULL

comp = rbind(cbind('H_28',  '-', 'C_28'),
             cbind('D_28',  '-', 'C_28'),
             cbind('HD_28',  '-', 'C_28'),
             cbind('W_28',  '-', 'C_28'))

for (i in 1:nrow(comp)) {
  
  # print(i)
  # print( as.character(comp[i, ]))
  
  dtX = data.long[data.long$mygroup %in% comp[i, c(1,3)] & data.long$Genotype %in% genotype, ]

  stress = dtX[dtX$mygroup == as.character(comp[i, 1]), ]
  Ctrl = dtX[dtX$mygroup == as.character(comp[i, 3]), ]
  
  StressC = gsub('_.*', '', comp[i, 1])
  CtrlC = gsub('_.*', '', comp[i, 3])
  
  log2FC = my.logFC(Ctrl = Ctrl, stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  p = my.customised.t.test(Ctrl = Ctrl, Stress = stress, StressC = StressC, CtrlC = CtrlC, genotype = genotype)
  
  De.log2FC = rbind(De.log2FC, log2FC)
  De.p = rbind(De.p, p)
  
  
}


Desiree.stat = merge(De.log2FC, De.p,
                     by.y = c("Genotype", "molecule", "group2", "group1"),
                     by.x = c("Genotype", "molecule", "stressGroup", "controlGroup"),
                     all = TRUE)


openxlsx::write.xlsx(Desiree.stat, '../reports/Desiree.logFC-metabolomics_tubers.xlsx', asTable = TRUE)

7.5 heatmap

comp = as.data.frame(comp)
forHeatmap = comp[intersect(grep('C', comp[, 1], invert = TRUE), grep('C', comp[, 3])), ]
forHeatmap$comparison = paste(forHeatmap[, 1], forHeatmap[, 3], sep = '-')

df = Desiree.stat[Desiree.stat$stressGroup %in% forHeatmap[, 1] & Desiree.stat$controlGroup %in% forHeatmap[, 3], ]

my.pheatmap.logFC(df, forHeatmap,
                  main = 'Desiree-tubers-metabolomics',
                  cluster_rows = FALSE, cluster_cols = FALSE,
                  subdir = 'tubers-metabolomics',
                  width = 6, height = 6, cnt = cnt)
# rm()
gc()
gc()
devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.4.1 (2024-06-14 ucrt)
##  os       Windows 10 x64 (build 19045)
##  system   x86_64, mingw32
##  ui       RTerm
##  language (EN)
##  collate  English_United Kingdom.utf8
##  ctype    English_United Kingdom.utf8
##  tz       Europe/Ljubljana
##  date     2024-11-05
##  pandoc   3.2 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package      * version  date (UTC) lib source
##  abind          1.4-8    2024-09-12 [1] CRAN (R 4.4.1)
##  assertthat     0.2.1    2019-03-21 [1] CRAN (R 4.4.1)
##  backports      1.5.0    2024-05-23 [1] CRAN (R 4.4.0)
##  broom          1.0.7    2024-09-26 [1] CRAN (R 4.4.1)
##  bslib          0.8.0    2024-07-29 [1] CRAN (R 4.4.1)
##  ca             0.71.1   2020-01-24 [1] CRAN (R 4.4.1)
##  cachem         1.1.0    2024-05-16 [1] CRAN (R 4.4.1)
##  car            3.1-3    2024-09-27 [1] CRAN (R 4.4.1)
##  carData        3.0-5    2022-01-06 [1] CRAN (R 4.4.1)
##  cli            3.6.3    2024-06-21 [1] CRAN (R 4.4.1)
##  codetools      0.2-20   2024-03-31 [2] CRAN (R 4.4.1)
##  colorspace     2.1-1    2024-07-26 [1] CRAN (R 4.4.1)
##  corrplot     * 0.95     2024-10-14 [1] CRAN (R 4.4.1)
##  crosstalk      1.2.1    2023-11-23 [1] CRAN (R 4.4.1)
##  data.table     1.16.2   2024-10-10 [1] CRAN (R 4.4.1)
##  dendextend     1.18.1   2024-10-13 [1] CRAN (R 4.4.1)
##  devtools       2.4.5    2022-10-11 [1] CRAN (R 4.4.1)
##  digest         0.6.37   2024-08-19 [1] CRAN (R 4.4.1)
##  dplyr          1.1.4    2023-11-17 [1] CRAN (R 4.4.1)
##  ellipsis       0.3.2    2021-04-29 [1] CRAN (R 4.4.1)
##  evaluate       1.0.1    2024-10-10 [1] CRAN (R 4.4.1)
##  fansi          1.0.6    2023-12-08 [1] CRAN (R 4.4.1)
##  farver         2.1.2    2024-05-13 [1] CRAN (R 4.4.1)
##  fastmap        1.2.0    2024-05-15 [1] CRAN (R 4.4.1)
##  foreach        1.5.2    2022-02-02 [1] CRAN (R 4.4.1)
##  Formula        1.2-5    2023-02-24 [1] CRAN (R 4.4.0)
##  fs             1.6.5    2024-10-30 [1] CRAN (R 4.4.1)
##  generics       0.1.3    2022-07-05 [1] CRAN (R 4.4.1)
##  ggh4x        * 0.2.8    2024-01-23 [1] CRAN (R 4.4.1)
##  ggplot2      * 3.5.1    2024-04-23 [1] CRAN (R 4.4.1)
##  glue           1.8.0    2024-09-30 [1] CRAN (R 4.4.1)
##  gridExtra      2.3      2017-09-09 [1] CRAN (R 4.4.1)
##  gtable         0.3.6    2024-10-25 [1] CRAN (R 4.4.1)
##  heatmaply    * 1.5.0    2023-10-06 [1] CRAN (R 4.4.1)
##  highr          0.11     2024-05-26 [1] CRAN (R 4.4.1)
##  htmltools      0.5.8.1  2024-04-04 [1] CRAN (R 4.4.1)
##  htmlwidgets  * 1.6.4    2023-12-06 [1] CRAN (R 4.4.1)
##  httpuv         1.6.15   2024-03-26 [1] CRAN (R 4.4.1)
##  httr           1.4.7    2023-08-15 [1] CRAN (R 4.4.1)
##  iterators      1.0.14   2022-02-05 [1] CRAN (R 4.4.1)
##  jquerylib      0.1.4    2021-04-26 [1] CRAN (R 4.4.1)
##  jsonlite       1.8.9    2024-09-20 [1] CRAN (R 4.4.1)
##  knitr          1.48     2024-07-07 [1] CRAN (R 4.4.1)
##  labeling       0.4.3    2023-08-29 [1] CRAN (R 4.4.0)
##  later          1.3.2    2023-12-06 [1] CRAN (R 4.4.1)
##  lattice        0.22-6   2024-03-20 [2] CRAN (R 4.4.1)
##  lazyeval       0.2.2    2019-03-15 [1] CRAN (R 4.4.1)
##  lifecycle      1.0.4    2023-11-07 [1] CRAN (R 4.4.1)
##  magrittr     * 2.0.3    2022-03-30 [1] CRAN (R 4.4.1)
##  memoise        2.0.1    2021-11-26 [1] CRAN (R 4.4.1)
##  mime           0.12     2021-09-28 [1] CRAN (R 4.4.0)
##  miniUI         0.1.1.1  2018-05-18 [1] CRAN (R 4.4.1)
##  mnormt         2.1.1    2022-09-26 [1] CRAN (R 4.4.0)
##  munsell        0.5.1    2024-04-01 [1] CRAN (R 4.4.1)
##  nlme           3.1-166  2024-08-14 [2] CRAN (R 4.4.1)
##  openxlsx       4.2.7.1  2024-09-20 [1] CRAN (R 4.4.1)
##  pheatmap       1.0.12   2019-01-04 [1] CRAN (R 4.4.1)
##  pillar         1.9.0    2023-03-22 [1] CRAN (R 4.4.1)
##  pkgbuild       1.4.5    2024-10-28 [1] CRAN (R 4.4.1)
##  pkgconfig      2.0.3    2019-09-22 [1] CRAN (R 4.4.1)
##  pkgload        1.4.0    2024-06-28 [1] CRAN (R 4.4.1)
##  plotly       * 4.10.4   2024-01-13 [1] CRAN (R 4.4.1)
##  plyr           1.8.9    2023-10-02 [1] CRAN (R 4.4.1)
##  profvis        0.4.0    2024-09-20 [1] CRAN (R 4.4.1)
##  promises       1.3.0    2024-04-05 [1] CRAN (R 4.4.1)
##  psych          2.4.6.26 2024-06-27 [1] CRAN (R 4.4.1)
##  purrr          1.0.2    2023-08-10 [1] CRAN (R 4.4.1)
##  R6             2.5.1    2021-08-19 [1] CRAN (R 4.4.1)
##  ragg           1.3.3    2024-09-11 [1] CRAN (R 4.4.1)
##  RColorBrewer * 1.1-3    2022-04-03 [1] CRAN (R 4.4.0)
##  Rcpp           1.0.13-1 2024-11-02 [1] CRAN (R 4.4.1)
##  registry       0.5-1    2019-03-05 [1] CRAN (R 4.4.0)
##  remotes        2.5.0    2024-03-17 [1] CRAN (R 4.4.1)
##  reshape2       1.4.4    2020-04-09 [1] CRAN (R 4.4.1)
##  rlang          1.1.4    2024-06-04 [1] CRAN (R 4.4.1)
##  rmarkdown      2.29     2024-11-04 [1] CRAN (R 4.4.1)
##  rstatix      * 0.7.2    2023-02-01 [1] CRAN (R 4.4.1)
##  rstudioapi     0.17.1   2024-10-22 [1] CRAN (R 4.4.1)
##  sass           0.4.9    2024-03-15 [1] CRAN (R 4.4.1)
##  scales       * 1.3.0    2023-11-28 [1] CRAN (R 4.4.1)
##  seriation      1.5.6    2024-08-19 [1] CRAN (R 4.4.1)
##  sessioninfo    1.2.2    2021-12-06 [1] CRAN (R 4.4.1)
##  shiny          1.9.1    2024-08-01 [1] CRAN (R 4.4.1)
##  stringi        1.8.4    2024-05-06 [1] CRAN (R 4.4.0)
##  stringr        1.5.1    2023-11-14 [1] CRAN (R 4.4.1)
##  svglite        2.1.3    2023-12-08 [1] CRAN (R 4.4.1)
##  systemfonts    1.1.0    2024-05-15 [1] CRAN (R 4.4.1)
##  textshaping    0.4.0    2024-05-24 [1] CRAN (R 4.4.1)
##  tibble         3.2.1    2023-03-20 [1] CRAN (R 4.4.1)
##  tidyr          1.3.1    2024-01-24 [1] CRAN (R 4.4.1)
##  tidyselect     1.2.1    2024-03-11 [1] CRAN (R 4.4.1)
##  TSP            1.2-4    2023-04-04 [1] CRAN (R 4.4.1)
##  urlchecker     1.0.1    2021-11-30 [1] CRAN (R 4.4.1)
##  usethis        3.0.0    2024-07-29 [1] CRAN (R 4.4.1)
##  utf8           1.2.4    2023-10-22 [1] CRAN (R 4.4.1)
##  vctrs          0.6.5    2023-12-01 [1] CRAN (R 4.4.1)
##  viridis      * 0.6.5    2024-01-29 [1] CRAN (R 4.4.1)
##  viridisLite  * 0.4.2    2023-05-02 [1] CRAN (R 4.4.1)
##  webshot        0.5.5    2023-06-26 [1] CRAN (R 4.4.1)
##  withr          3.0.2    2024-10-28 [1] CRAN (R 4.4.1)
##  xfun           0.49     2024-10-31 [1] CRAN (R 4.4.1)
##  xtable         1.8-4    2019-04-21 [1] CRAN (R 4.4.1)
##  yaml           2.3.10   2024-07-26 [1] CRAN (R 4.4.1)
##  zip            2.3.1    2024-01-27 [1] CRAN (R 4.4.1)
## 
##  [1] C:/Users/majaz/AppData/Local/R/win-library/4.4
##  [2] C:/Program Files/R/R-4.4.1/library
## 
## ──────────────────────────────────────────────────────────────────────────────